IRPHE/1997-LRV 

Hamiltonian Dynamics and the Phase Transition of the XY Model 



Xavier Leoncini, and Alberto D. Verga,* 
Institut de Recherche sur les Phenomenes Hors Equilibre,^ 
12, avenue General Leclerc, F-13003 Marseille, France. 



and Stefano Ruffo 
Dipartimento di Energetica "S. Stecco", 
Universitd degli Studi di Firenze, INFN and INFM, 
Via di S.Marta, 3 50139 Firenze, Italy. 
(February 1, 2008) 



Abstract 

A Hamiltonian dynamics is defined for the XY model by adding a kinetic 
energy term. Thermodynamical properties (total energy, magnetization, vor- 
ticity) derived from microcanonical simulations of this model are found to be 
in agreement with canonical Monte-Carlo results in the explored temperature 
region. The behavior of the magnetization and the energy as functions of 
the temperature are thoroughly investigated, taking into account finite size 
effects. By representing the spin field as a superposition of random phased 
waves, we derive a nonlinear dispersion relation whose solutions allow the 
computation of thermodynamical quantities, which agree quantitatively with 
those obtained in numerical experiments, up to temperatures close to the 
transition. At low temperatures the propagation of phonons is the dominant 
phenomenon, while above the phase transition the system splits into ordered 
domains separated by interfaces populated by topological defects. In the high 
temperature phase, spins rotate, and an analogy with an Ising-like system can 
be established, leading to a theoretical prediction of the critical temperature 
Tkt « 0.855. 
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I. INTRODUCTION 



The two dimensional XY model, also known as the planar spin model, presents many 
interesting behaviors. Despite the presence of a continuous symmetry group, a particular 
form of phase transition exists , which can be characterized by the change in the behavior 
of the correlation functions. In the low temperature phase these latter have power law decay, 
showing that the system is in a long range order state; while they decay exponentially at high 
temperatures, the long range order is broken, even though thermodynamic quantities remain 
smooth across the transition 0. These observations have been interpreted by Kosterlitz and 
Thouless [Q using an analogy with the transition of a Coulomb gas from a dielectric phase, 
where charges are bounded into dipoles, to a plasma (conducting) phase where temperature 
fluctuations destroy the dipoles, and the charges become free. In the XY system the charges 
are replaced by topological excitations called vortices. 

From the analytical point of view, beyond the spin- wave approximation [^, and the 
Villain model 0, the use of renormalization group in the critical region has been the main 
issue In order to confirm the analytical results and to satisfy the need for a better 

understanding of the transition, many numerical studies have been performed. Different 



dynamics, like Monte-Carlo [y] or Langevin ||l0[, have been introduced. The system has been 



confirmed to be dominated by spin wave excitations in the low temperatures region. The 
transition region is then determined using the information on the correlation functions, and is 
interpreted in terms of dipoles unbinding. However, there is a number of observations in the 
literature that seems to complicate this simple picture of the phase transition mechanism. 
Among the different results those worth mentioning are: the visualization of the vortex 
distribution 0, the presence of domains delimited by topological defects [0, the precise 
determination of the transition temperature T^t ~ 0.89 much lower than the transition 
temperature of the Villain model. Results are also found on the interaction between vortices 



p!0| , or the vortex interaction energy ||T2 |. 



Although the basic mechanism of the Kosterlitz-Thouless transition, in terms of the 
breaking of vortex dipoles associated to the emergence of a disordered state, is well under- 
stood, the observation of the spatial distribution of defects, which is not uniform (defects 
tend to appear organized into clusters at temperatures slighty larger than the transition tem- 
perature), and the presence of large ordered domains where the spins are almost parallel, 
seem to indicate that the physics of the phase transition is not exhausted by this unbinding 
process but that some kind of partial local order is present even beyond the transition tem- 
perature. The investigation of the system properties near the transition is one of the points 
addressed in this paper. 

On the other hand, the question whether statistical physics is able to describe over 
a wide range of temperatures the behavior of a classical Hamiltonian system with many 
degrees of freedom, still remains an open one (see e.g. the literature on the Fermi-Pasta- 



Ulam model quoted in [|T3[)- Therefore, it seemed interesting to consider the XY model from 
the point of view of Hamiltonian dynamical systems by adding a kinetic energy term to the 
XY Hamiltonian. Such kind of approach has proved to be fruitful of interesting information 



on typical relaxation time scales and collective behaviors in the ID case |I4| and in the 



mean field approximation |T^. More recently Lyapunov exponents have been computed. 



confirming the presence of long relaxation times to equilibrium both in the very low and 
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in the very high temperature hmits in one dimension [|T^] and the presence of a change in 
slope of the maximal Lyapunov exponent vs the energy near the transition temperature in 
2D (see also Ref. [|I^ for a preliminary study of this latter phenomenon). 

The present model is one of coupled rotators (spins) sitting on the sites of a square 
lattice interacting with near neighbors, whose statistical properties are described by the mi- 
crocanonical ensemble, the total energy being set by the initial conditions. Convergence to 
a Gibbsian equilibrium distribution is not given for granted in the whole temperature range 
on accessible time scales. Therefore, a particular attention must be devoted to the temporal 
behavior of the different quantities which characterize macroscopic thermodynamical prop- 
erties. One notices, however, that the spatial topology is unchanged and that the existence 
of a continuous symmetry group of rotation is also maintained in the dynamics. This implies 
that a Kosterlitz-Thouless type transition must be observed (in the thermodynamic limit) 
with typically strong finite size effects, like the existence of a non-zero magnetization []5|,p!9|. 

In this paper we concentrate on a study of those properties of the dynamics of the XY 
model which reproduce equilibrium features, postponing to a future work the study of non 
equilibrium effects. We anticipate that we are able to reproduce most of the equilibrium 
behaviors of macroscopic quantities, which makes microcanonical dynamics competitive with 
canonical Monte Carlo, since we have not here to extract random numbers, stochasticity 
being supplied by the intrinsic chaoticity of the model. 

The fact that we actually deal with a dynamical system, allows us to develop an original 
analytical approach to the study of the thermodynamics, which is based on the approximate 
solution of the equations of motion. This method is based upon the ergodic properties of 
the dynamics and on the separation of temporal scales present in the spin motion. 

In Section |I| we present the Hamiltonian model, the basic aspects of the numerical com- 
putations and we introduce the thermodynamic and dynamical quantities that characterize 
the state of the system. In Section |T| we compare the statistical properties of our model 
to the usual ones, using both analytical and numerical approaches. We derive the ther- 
modynamical properties using the hypothesis that the dynamics is dominated by random 
phased waves. We then propose, in Section |^ an Ising-like model based on the observation 
that, above the transition, synchronized regions of spins appear. The relation of the XY 
with an Ising-like model allows us to describe the high temperature phase and to derive an 
approximate value of the critical temperature. A brief summary of the main results and 
conclusions is presented in Section 0. 



II. THE HAMILTONIAN MODEL, BASIC PROPERTIES AND NUMERICAL 

COMPUTATIONS 

The XY model was introduced in statistical mechanics as a two-dimensional version of 
the Heisenberg Hamiltonian. The spins are fixed on the sites of a square lattice, and are 
characterized by a rotation angle 9i G [— tt, tt] 

N 

HxY = JY.\^-cos{ei-e,)] , (2.1) 

where J is the coupling constant (with J > corresponding to the ferromagnetic case, that 
we study here), i and j label the sites of a square lattice of side \fN , i.e. i = (ix^iy) 
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with 1 < ix,iy < \/N of coordinates {x,y). The summation is extended over all i and its 
neighboring sites j. In the following, without loss of generality, we set J = 1 and the lattice 
step to unity. 

The spins evolve in time, 6i = Oi{t), after adding a kinetic energy term to the XY 
Hamiltonian, 

H = J:^ + Hxy, (2.2) 
i=i ^ 

where = Qi is the spin momentum. The choice J = 1 is equivalent to set time units 
and to rescale accordingly momentum (in these units the "inertia" is also unity). With this 
kinetic energy term the spins become in fact rotators, and the XY model becomes a system 
of coupled rotators. The equations of motion are 

4 

e^{t) = -Y.^m[elt)-e,{t)\^ (2-3) 

where the summation is over the four neighbors j of site i. In addition to the energy H = E, 
there exists a second constant of the motion, the total angular momentum P = J^iPi^ which 
can be chosen to be zero. We choose periodic boundary conditions in both the x and the y 
directions. Numerical integration of Eq. (|2.3| ) is performed using the Verlet algorithm, which 
conserves the energy 0{At'^), At being the time step, but exactly preserves momentum and 
the symplectic structure. 

Thermo dynamical quantities are computed by averaging over time and over the sites 
a single orbit (the evolution of the system from a given initial condition). Typically, the 
system is started with a Gaussian distribution of momenta and with all the spins pointing 
in the same direction 6i = 60, for reasons that are clarified in the following. No strong 
dependence on the chosen initial condition in this class was observed, but one could improve 
statistically our results by averaging over many orbits with different initial conditions and 
the same energy. 

The temperature is computed through the average squared momentum per spin 

1 ^ 

T=^T.P^ity^ (2-4) 
i=i 

where the overbar stays for temporal averaging. 

The thermodynamical state can be characterized by several macroscopic variables: the 
internal energy per spin h = h{T) = E{T, N)/N {E being the constant total energy of 
the system); the magnetization M = M(T, A^), which, as mentioned in the Introduction, 
vanishes for — 00, but is sizeable for any finite A^; the density of topological defects py, 
or vortices, which is intimately related to the mechanism underlying the phase transition. 

The magnetization M = (Mj., My) is given by 

M = M{t) = ^ E (cos Oi, sin Oi) , (2.5) 



4 



and describe the mean orientation of the spin field. A topological defect is identified, as 
usual, by computing the total angle circulation on a given plaquette (the sum of the four 
plaquette angles mod[— vr, vr]); when it equals ±27r this quantity identifies a positive (nega- 
tive) unitary vortex on the plaquette. The total density of vortices, or vorticity, which is an 
intensive quantity, is then given by 

^« = ^ - ^j) mod [-vr, tt] , (2.6) 

where denotes the sites j of the i-th plaquette. We have studied the temporal evolution 
and the spatial distribution of the vortices as a function of the temperature. Because of the 
periodic boundary conditions and since P = 0, the number of positive vortices equals the 
number of negative ones. 

We have performed simulations using various sets of parameters, to study the system 
behavior depending on the temperature, the number of spins, and the total time to test the 
stationarity of the relaxed state. Having chosen parallel spins, one has an initially vanishing 
potential energy, allowing the exploration of the low temperature region through the reduc- 
tion of kinetic energy. Random angles in the [— vr, vr] interval would give an energy per spin 
of 2, which would then remain fixed producing a high temperature configuration, as can be 
also deduced from the h{T) function we compute below (see Fig. |l]). Therefore, initially, 
the magnetization is |M| = 1. Other initial conditions, with |M| =0 were investigated, 
for instance with half the spins oriented in the x direction, and the other half in the —x 
direction. After a transient, the system relaxes to the same thermodynamical state. 




0.5 - # 
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FIG. 1. Total energy per spin as a function of the temperature. At low tempera- 
ture the energy grows as T and at high temperature the energy tends to T/2 + 2. The 
Kosterlitz-Thouless phase transition occurs at Tkt ~ 0.89. +, x, and * signs refer to the 
Hamiltonian dynamical simulations for different lattice sizes, = 64?, N = 128^, and 

= 256^ respectively; circles refer to canonical Monte-Carlo simulations for a lattice size 
of N = 100^. 

In Fig. |l|, we plot the total energy per spin h{T). The low (T ^ 1, T 1 corresponds to 
the value for which the kinetic energy is of the same order of the potential energy) and high 
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(T ^ 1) temperature behaviors are easy to understand. At low temperature equipartition 
of kinetic and potential energies gives h{T) ?a T (this is the linear regime, where angles 
between neighboring spins are small). At high temperature, angles are uniformly distributed 
in [— '7r,7r], the cosine interaction in (p.2|) is negligible with respect to the kinetic energy, 
then h{T) ^ T/2 + 2. These simple arguments suggest that for T <^ 1, the spin field 
is almost linear and can be represented as a superposition of waves, which correspond to 
phonons, while for T 3> 1, the potential energy is negligeable, and this field becomes a set 
of almost free fast rotators. The change in the behavior of the h{T) curve starts around the 
Kosterlitz-Thouless critical temperature T^t- The peak in the specific heat, related to the 
second derivative of the energy, occurs instead at a somewhat higher temperature, a fact 
well documented in Monte-Carlo computations. 



* 
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FIG. 2. Absolute value of the magnetization as a function of the temperature. Finite 
size effects give a finite magnetization at low temperatures. At low temperatures the mag- 
netization decreases only logarithmically with the number of spins. -|-, x, and * signs refer 
to the Hamiltonian dynamical simulations for different lattice sizes, 

and N = 256^ respectively; circles refer to Monte-Carlo simulations for a lattice size of 
N = 100^. 



Figure ^ shows the absolute value of the magnetization as measured using formula (p. 51). 
Although in the thermodynamic limit the magnetization must vanish, in a finite system and 
for low temperature we expect an observable macroscopic magnetization, which decreases 
logarithmically with the number of spins as noted by Berezinskii [Q. This slow decrease 
with system size is hardly observable in Fig. ^ but points corresponding to larger sizes give 
systematically a smaller magnetization. In the high temperature region, since angles are 
randomly distributed in space at any time, the magnetization vanishes algebraically with 
the number of spins. 

The behavior of vorticity with temperature, plotted in Fig. ^, is related to the form 
of the h{T) curve. Indeed, at low temperature aligned angle configurations are typical and 
consequently vorticity vanishes, p^, —>■ 0. In the high temperature regime, spins are randomly 
distributed and the vortex density approaches the asymptotic value — >■ 1/3. A dramatic 
growth of the vorticity occurs in the phase transition region (the temperature in the range 
T ^ 0.8 to 1.5). 



6 



FIG. 3. Number density of vortices as a function of the temperature normalized to 
the total number of spins. This quantity is intensive. +, x, and * signs refer to the 



Hamiltonian dynamical simulations for different lattice sizes, 
N = 256^ respectively. 



642, iV = ^2S^ g_^j 



All our results are in perfect agreement with those of Monte-Carlo simulations (see e.g. 
), some of them are plotted in the Figs. and |^ for comparison. Therefore, microcanonical 
(Hamiltonian dynamics) and canonical (Monte-Carlo) computations give, at least for the 
class of initial conditions studied here, the same thermodynamical equilibrium states. 



III. DYNAMICAL AND STATISTICAL PROPERTIES 

We know that at low temperatures, the main contribution to the canonical partition 
function comes from the configurations in which the spins are almost aligned, configurations 
where the angle differences 6i — 6j are small. In such a situation, the equations of motion 
( p.3|) can be linearized, and the spin field is therefore represented as a superposition of 
linear waves. We will here determine the effective dispersion relation for these waves at 
low temperature using consistency relations with temperature and internal energy. Let us 
introduce a representation of the spin angles in the form of a random phased field 

ei = J2'^k cos(^^) , i^l = k ■ Xi - Ukt + (f)k , (3.1) 

k 

where the summation is over the wave-vectors k = (kx, ky) = 2Tc(nx, Uy) / V^, with n^., Uy = 
1, . . . , a//V integers; the wave spectrum is given by the a^, and the phases (pk are supposed 
to be random, uniformly distributed in the circle. We also should have added to (|3.1| ) a term 
of the form reflecting the individual rotation of spins with some frequency fij, which in 
fact can be considered a function of the temperature. However, at low temperature most of 
these frequencies must vanish, because the energy needed to trigger the free rotation is only 
reached at energies of the order of h = 2 (or temperatures of the order of T = 1.5, as can 
be seen in Fig. |I|), when a spin, considered as a perturbed pendulum, crosses the separatrix. 
Then, in the following we put fij = in the low temperature (this term becomes important 
at high temperatures). 
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Although the main assumption in ( |3.1| ), that is the random character of the wave phases, 
may only be justified a posteriori, we performed some numerical tests which showed that 
this ansatz is consistent with the properties of the system (at low and also at finite temper- 
atures). We found that the probability distribution of the spin velocities is almost uniform, 
and that moreover, the motion of one single spin is ergodic (in the sense that its temporal 
mean and variance coincide with the spatial ones). In addition, if the system were linear this 
assumption is equivalent to the assumption of thermodynamic equilibrium (with a bath at 
temperature T). In fact this is the case at very low and very high T; of course, in the inter- 
mediate range the infiuence of the nolinearity becomes important, for instance, establishing 
non-local interactions or some type of self-organization, in which case this assumption would 
not be necessarly valid. 

In formula ( |3.1| ) there are two unknown functions of the wave-number which remain to 
be determined, the wave frequencies Uk and the spectrum a^. Note, in addition, that ( p.l|) 
is in fact a change of variables of the angles 6*, to the amplitudes ak, where implicitly, a slow 
temporal dependence may be included, the fast time dependence being assured by the phase 
ujkt. In this sense, formula ( |3.1| ) is rather general and can take into account a large variety 
of field states. 

The definition of the temperature imposes a constraint on the form of the wave spectrum. 
Indeed, replacing formula ( |3.1|) into the definition of the temperature (|2.4| ) we obtain 

1 ^ • 

i 

= ^E«^^L (3-2) 
^ k 

where (■ ■ ■) stays for averaging over the random phases (j)k-, and where we used the identity 

(cos(V^^) cos(V^^,)) = \^k,k' , 

being 6k,k' the Kronecker symbol. For instance, if equipartition of the energy among different 
degrees of freedom is assumed, one obtains the usual Jeans spectrum, given by 



2T 



= ^ , (3.3) 

which means that each degree of freedom takes a fraction T/2 of the total kinetic energy. 
In fact, as we will show below, the detailed knowledge of the spectrum turns out not to be 
necessary to derive the dispersion relation and compute thermodynamical quantities. 

Before proceeding to the computation of the dispersion relation, let us derive a general 
expression of the energy per spin h, 

rp 1 N 

h = - + 2--Y. {oos{9. - %)) , (3.4) 

<i,j> 

where the first term comes from the definition of the temperature, the second term is the 
constant added for convenience to the Hamiltonian to make the energy vanish at zero tem- 
perature, and the last term includes, for each lattice site, a summation over two neighbors. 
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for instance west (in the x direction) and south (in the y direction). When the expansion 
( p.l|) is introduced into Eq. (|3.4| ) the typical contribution to the sum in the third term on 
the r.h.s. is 



U 



{id) 
N 



COS 



N 



.n=l 



where n relabel the modes and 



A 



B 



2 sin 



1X7, X 'i 



(ij) k ■ (Xj + Xj) 



are local in n and independent on the phases (pkn- Let us split the summation over the wave- 
numbers into two terms, one containing the phase, 0a:jv) ^"^^ the other with the remaining 
n = 1, . . . , — 1 phases. Since 



cos 
sin 



4f ^HBt'^' + 

(id) 

'N 



4f MBi:f + 



N 
JV 



(3.5) 



where Jn is the Bessel function of zero order, we obtain the recursion relation 



U 



N 



which readily gives the result 
We finally get 



U 



N 



^ k k 



(3.6) 



where (5^ = 2afcSinfcj,./2 and jSy = 2akSinky/2 and the last two products come from 
north/south and east/west neighbours on the lattice. Eq. ( |3.6| ) is the general expression 
for the energy per spin, as long as the system is dominated by random phased waves. We 
note that in the case of a symmetric spectrum ak^^ky = Ciky,k^ (this is the case for the Jeans 
spectrum, when the dispersion relation satisfies u!k^,ky = ^ky,k^), the products of zero order 
Bessel functions are equal. 



7, = I[Mf3.)=llMPy), 

^ k k 



(3.7) 



where v is related to h by h{T) = T/2 + 2 — v(T), and is the average potential energy. In 
order to obtain an explicit form of the energy per spin as a function of the temperature, 
h = h{T), we must determine in a self-consistent way the dispersion relation; this may 
be done by using the equation of motion in its full nonlinear form. But let us begin with 
the linear case, which will serve as a guide to the self-consistent computation. The linear 
equation of motion reads 
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Substituting one obtains 

-J2 ^lakcos^l 



- 

<j> 



-E/5.[sin(^^ + y)-sin(^^-|)] + 

k 
k 

-4^«.[sin^(y) + sin^4)] cos^^ 



(3.8) 



(3.9) 



although the hnearity of the equation allows an identification term by term in the summation 
over k, a different strategy for the solution would be to average (|3.9| ) over all the phases 
but one, e.g. over the phases 0^/ ^ (pk after multiplying both sides by cos?/'^, in order to 
isolate only one term in the summation on the left hand side. One finally obtains the linear 
dispersion relation 



l^fc — ^Ok 



h 

4(sm h sm 



2 ^y- 



(3.10) 



The frequency spectrum is symmetric with respect to the exchange ^ ky] in the following 
we assume, without loss of generality, that the spectrum is symmetric with respect to this 
transformation. 

This procedure can be generalized to the nonlinear case, although the r.h.s. of Eq. ( [^.3[ ) 
is no more local in 0^. After substituting the expansion ( p.l|) into Eq. (|2.3| ), we must average 
over the phases (pk' 7^ 4>k, which we denote <■■■>', terms of the form 



{sm{e - Oe))' = - ^sin 



k 

/3xSin(V^fc + y) 

. k 



(3.11) 



where 9e is e.g. the east neighbor of angle 9, having a phase ipE = i^k+k^ = k-x+k^—Ukt+cpk- 
We now split the summation in ( p.ll|) into a term containing the phase 0^ and the others, 
and develop the sine of the summation of two terms into a product. The averaging of these 
two terms reduces to 



cos 



E/?.sin(^;±f) 

- k' 



(3.12) 



because the averaging of the sine term is zero. The two neighbors in the x direction (west 
and east sites) give the following expression 



(sin(^^ - Be))' + {sm{e - Ow))' 



sm 



(3^ sin(V'fc + y ) 



sm 



PxSiniipk - y) 



(3.13) 



Those in the y direction give exactly the same contribution with (3^ (3y. On the other 
hand, the left hand side of Eq. ( p.3|) , the temporal part of the equation of motion, is reduced, 
after substitution of the random waves expression, to 
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(^Oi'j = -ulak COS ipl 



(3.14) 



As before, in the linear computation, we multiply both sides of the equation of motion ( p.l4| ) 



and ( 3.13| ) (including the similar terms for the north and south neighbors), by cos('?/^^) and 



average over the phase 0^, noting that this average involves first order Bessel functions Ji. 
We finally find the desired dispersion relation 



ujlak = 2v 



smi — ) — -— — + smi ^ 



2 ' Jo(/5,) ' 2 ' Jo{(3y 



(3.15) 



This dispersion relation is nonlinear, i.e. the frequency depends on the spectrum amplitudes, 
in two ways: implicitly through the a^s in v and explicitly in the Bessel functions. To go 
further, let us investigate the dependence of the arguments of the Bessel functions on the 
parameters T and A^, and the related limiting form of the total energy per spin h. We 
know that v, being an intensive thermodynamic quantity, does not depend on the number of 
spins (the size of the system); moreover, the linear frequency (|3.1CI| ) is bounded from below, 
ciJofc ~ siia^{kx/2) > 0{1/N). We also note, as it is natural for a system near an equilibrium 
state, that a large number of k modes must contribute to the energy of the system, and 
then in general we have alsmkx/2 —>■ when N oo, for a large range of temperatures. 
If this were not the case, the energy would be concentrated in a few high k modes, which 
is clearly in contrast with the observations. For the Jeans spectrum, one finds specifically 
alsmkx/2 0{T / \/N). Using these approximations we can now develop the logarithm of 
the product of the Bessel functions in (|3.7| ) to obtain 



t; = 2exp|--^Q;^^^^| . (3.16) 

In the same approximation, for fixed temperatures and large lattices, the dispersion relation 
( p.l5| ) reduces to 



^1="-^- (3.17) 

Introducing this expression into the formula ( p.l6| ) for v and using the definition of the 
temperature in terms of random phased waves ( p.2| ) we obtain an implicit equation for the 
potential energy. 



V = 2 exp 



T 
2^ 



(3.18) 



We notice that the nonlinear dispersion relation reduces to the linear one for T = 0. The 
nonlinear effects, taken into account in (|3.17|) , appear as a renormalization of the phonon 
frequency (energy) due to the coupling of the phonons with a thermal bath created by the 
other phonons like in a mean field. We also note that the factor 2 found in the right hand 
side of ( |3.18| ) and in the exponent can be related to the lattice dimensionality (it results 
from the addition of the x and y terms). Moreover, as we anticipated, formula ( |3.18|) does 
not depend explicitly on the form of the spectrum a^. 
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Temperature 



FIG. 4. Plot of the energy density versus temperature computed from Eq. ( p. 18 ). The 
straight hne characterize the random field state, the lower branch, the phonons (spin 
waves), and the upper branch, lattice excitations of the vortex dipole type. 



Let us now investigate the implicit equation (|3.18|) by solving it numerically. In Fig. ^, 
where the energy density h versus the temperature T is plotted using formulas h = T/2 + 
2 — v{T) and (|3.18|) , we see three branches according to the solutions of ( 3. IS ). Although the 
three branches are in principle acceptable solutions, they would have different "weights" . For 
instance, if one associates to each branch a Gibbsian probability, proportional to exp(— /ij/T), 
where hi denotes the branch energy i = 0, 1, 2, the lower branch has, at a given temperature, 
the higher probability. The choice of a Gibbsian probablity distribution of states is not 
in contradiction with our microcanonical approach. Indeed, the energies hi are related to 
macroscopic equilibrium states, and the Gibbs probabilities represent the best choice, taking 
into account the usual constraints of maximal entropy and normalization. 

The straight line corresponds to the solution v = 0, which represents a random spin 
field. Indeed, in the high temperature region, since the potential energy is bounded, we can 
expect, that the rotators are freely rotating without any order, i.e. that {6i — 6j) is random 
with constant distribution in [0,27r). This gives us an energy density h = T/2 + 2 and 
v{T) = 0. The implicit solution of ( p. 18 ) extrapolates thus the good asymptotic behavior 
of the energy for high temperatures to the low temperature region. This solution exists 
for any temperature, but in the low temperature region spin configurations associated to 
this branch should be unstable (in the sense that for general initial conditions, as in our 
numerical computations, the system cannot evolve to this state). 

Below a certain temperature T = we also see that two other branches appear. As the 
lower branch has the higher probability, in the low temperature region, we may expect the 
lower energy branch, which we call the phonon branch, to be the physical relevant solution, 
and the third branch to have no physical meaning at all. However, if we now consider that in 
fact the lattice is populated with two types of species of different energies hi{T) and h2{T) 
corresponding to the lower (phonons) and upper (which we call vortex) branch respectively, 
we obtain 
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h 



(3.19) 



g-2/ii/T _^ g-2/12/T ' 

which is almost in total agreement with the numerical data up to T^,, as we can see in Fig. |^. 




0.5 1 1.5 

Temperature 

FIG. 5. Plot of the energy density versus the temperature. The circles refer to the 
numerical data, the dashed line to the phonon branch of Eq. ( 3.18| ), and the solid line to 
the theoretical averaged density of energy, Eq. ( p. 19 ). 



The upper branch can be related to "dipole vortices" since it starts at an energy equal to 
2 (see Fig. H) corresponding to a lattice composed by topological defects forming a perfect 
crystal {Oi — 6j = ±7r/2 and sin(6'j — 9j) = overall the lattice). This crystal state, 
which is in fact a periodic array of dipoles, is a stationary solution of the equations of motion 
at T = 0, around which no linear development can be made (the linear term is zero, vortices 
can be considered as nonlinear particles). The factor 2 in the exponents of ( p.l9| ) takes into 
account that vortices come in pairs (dipoles at low temperature) on the lattice in order to 
conserve the total circulation. 

The temperature determines the upper border of the wave dominated regime, above 
this temperature the wave branch disappears (as well as the vortex dipole branch). In order 
to determine precisely it is useful to rewrite (|3.18|) , as a formula for the temperature in 
terms of v, 



T = -2i;log(i;/2) 



(3.20) 



If we now calculate the extremum of T{v) for < f < 2, we obtain for v = 2/e = T^/2, 
T = = 4/e ~ 1.47, which is exactly the critical temperature found in using a 



Hartree-Fock approximation. This temperature also appears with the use of renormalization 
techniques in p2[. In this case, we note that the equation for v is formally identical to the one 



for the renormalization constant in this renormalization group calculation. However, the two 
methods are basically different, and in particular the relation between v and the dispersion 
relation (|3.17|) , allows us to get more precise informations on the different thermodynamical 
quantities of the system. For instance, as mentioned above, the relation between v and h, 
which is different from the relation between the renormalization constant and the energy, 
leads to a very good description of the h{T) curve up to T^. 
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FIG. 6. Finite size effects on ttie magnetization. Tlie upper plot sliows 
— 47rt;log(M)/log(2A^) as a function of the temperature; the hne refers to formula ( p.23| ); 
+, X, and o refer to numerical data for lattice sizes of N = 64?, N = 128^, and = 256^, 
respectively. The lower plot shows M{T), circles are for a N = 256^ lattice; the thin line is 
the analytical result and the dotted line is the result using only the contribution of linear 
spin waves: M = exp(— nog(2A'')/87r). 



Let us now evaluate the magnetization for a finite size system with periodic boundaries 
using the same approach as before: we assume that the spin field can be represented by a 
superposition of random phased waves We call 9o, the average over the lattice of the 

6i, which is a constant of motion, since the total momentum is conserved. An average of the 
magnetization over the random phases brings up to a similar calculation as the one for the 
energy per spin and leads to, 

(M) = n^oK)(coseo,sin^^o) (3.21) 

k 

This expression remains exact as long as the random phases approximation is valid. In 
order to develop the logarithm of this expression, we have to take into account that < 
0{T), which implies that the development is only valid in the low temperature regimes T -C 
1. Moreover a detailed knowledge of the spectrum is also required. Therefore, considering 
the observed almost flat spectrum of the momentum, we assume the equipartition of the 
kinetic energy among the modes and use the Jeans spectrum given in ( |3.3| ). The absolute 
value of the magnetization, (|M|), is given by the product over the different k of the Bessel 
functions. In the low temperature regime using ( p.l7| ) and (|3.3|) , we obtain for its logarithm, 

log ((IMI)) = - y ^ = y^ = --G (0) , (3.22) 

where G is the Green function of the linear wave equation, G{r) = I]fc6xp(iA; ■ r)/u!Qf^, with 
G{0) = (l/47r) log(2A^). The expression for the magnetization is then 

1 4^ 

(M) = (^— ) (cos ^^0, sin ^o) . (3.23) 

A plot of this expression, scaled in order to obtain a function of the temperature, 
— 47rf log(M)/ log(2A^), is shown in Fig. |^ (top). Within the errors of the numerical data. 
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the points collapse, at low temperature, to a unique curve. The result ( p.23|) , obtained 



including the correction due to the nonlinear contribution to the dispersion relation, sub- 
stantially improves the usual estimation based on the linear wave approximation (see the 
bottom of the figure). The agreement of the theoretical results and the numerical data can 
be made more precise by taking into account the different energy branches. Indeed, if we 
now take into account that the lattice is populated by two types of species, we can consider 
that only the phonons contribute to the magnetization, vortex dipoles having total magne- 
tization zero. As the density of phonons is given by rii = e-2'i-i/T^^g-2ft,i/T _|_ ^-2h2/T-^^ ^^^e 

observed magnetization must then be simply tt-iM. The absolute value of this quantity is 
plotted in Fig. ^ (bottom). The agreement with numerical data is valid up to the critical 
temperature. 

In the same way we compute the averaged absolute value of the magnetization 

(|Mp) = ^E(cos(e^-%)), 



and we obtain the expression 

i,j k 



(|Mr) = ^En^o(2a.sin(fc.^^)) , (3.24) 



which, in the low temperature regime, assuming a Jeans spectrum, and making the same 
expansion as in (p.22|) , leads to 



iMR = ]iii:-p(-^G 



where = Xi — Xj. A first order calculation leads to the same result as ( |3.23| ), implying 
that the variance Var(M) = (|Mp) — |(M)p vanishes at leading order in the temperature. 
However, at finite temperature, the two values can differ, and the variance of the magne- 
tization can have a nonzero value. This variance can moreover be used to characterize the 
phase transition, since it is the average over the lattice of the correlation function (e^'~^^ ); 
it is also related to the susceptibility x = (iV/T)Var(M). Indeed, in Fig. |^ we plot the 
magnetization variance as a function of the temperature, and find that up to temperatures 
of the order of the Kosterlitz-Thouless T^t, the variance remains very small, and suddenly 
it grows around T^t- This behavior suggests that the magnetization is distributed almost 
like a delta function, in the low temperatures regime (all the spins are pointing to the same 
direction), and that near the transition, the amount of randomly oriented spins increases 
dramatically. Higher order effects on the variance of the magnetization can be computed 
using the identity 

oo 

Jo{2zsma) = Jq{z) + ^ J^(z) cos(2mQ;) , 

m=l 



which gives, after replacing it into (|3.24|) 



Var(M) = |(M)|2 | ^ EH [1 + Mk{nj)] - 1 

i,j k 
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where 



Mkij-ij) = '^^^^''l cos {mk ■ Tij^ 



To the first nonvanishing order (m = 1) we obtain 



Va.(M)=2Q l(M)|^i:4 



This allows us to compute an approximate value of the suceptibihty 

T 

A similar result can be found in |E3i, in the limit v ^ 2. 



OOO-O-O- O- -© -G- - 



Temperature 

FIG. 7. Plot of the variance of the magnetization versus the temperature (A^ 
The transition region (T ~ 1) is clearly identified. 



128^ 



Concerning the density of vortices, we can see in Fig. ^ that around T ^ 1.3 the curve 
has an inflexion point. This change in the number of defects may be explained as follows. 
Let us consider a single plaquette, a defect may appear when the four angles are in increasing 
(decreasing) order, and the last angle is larger (smaller) than tt (— tt), if the first angle is 
set to zero. In our model, due to the continuous symmetry group, no particular direction is 
favored, and therefore the successive differences between the angles 6i—6j (discrete gradients) 
are all equivalent. A defect is then obtained when these gradients satisfy (in average) 
{\0i — 9j\) > 7r/3. At low temperatures the amplitudes of the local gradients are determined 
by the phonons and become steeper at higher temperatures, reaching a point where they 
are large enough to generate the vortices. Using the previous results, we calculate the 
temperature for which this condition is reached. 



2 ^ sin^ 



knr 



T 



v{T) 



71 

~9 



(3.25) 



If we substitute now this result into equation ( p. 181) , we find 
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T=— exp(^--j^l.27, (3.26) 

which is in very good agreement with the position of the inflexion point in the numerical 
data. One may consider that around this temperature a prohferation of vortices should 
occur. 



IV. THE PHASE TRANSITION 

The usual physical picture of the Kosterlitz-Thouless transition is based on the unbind- 
ing of vortex pairs: below the critical temperature long range correlations are established 
by spin waves, the phonons in our Hamiltonian model; above the critical temperature long 
range order is destroyed by the proliferation of free vortices. This simple physical picture of 
the mechanism of the transition does not exhaust the complex processes related to the ap- 
pearence of a kind of self-organization in the system: the vortex distribution is not uniform, 
vortices form clusters which separate domains of relatively well ordered spins. In this Section 
we study, using numerical simulations, the spatial distribution of topological defects and the 
equilibrium configurations of the system, using an Ising-like model, in order to describe the 
partial order which is still present above T^t- 

In the previous Section we showed that the ordered state of the system is dominated 
by phonons, and that a second type of excitations, the vortex dipoles, is also present. The 
energies of the two branches meet at a temperature Ty^, above which both modes disappear, 
suggesting that a phase transition to a disordered state, might occur. Within this picture, 

would be thus the temperature where the unbinding of vortex dipoles appears, destroying 
the long range order created by the waves. However, the actual transition temperature Tkt 
is much lower than T^, indicating that other processes, not included in the representation 
( ^TTD of the spin motion, take place. 

Indeed, in this low temperature analysis we have omitted the rotation of the spin, exclud- 
ing terms of the form flit. Moreover, most significantly, we have assumed that the phases of 
the spin field were random and uncorrelated from site to site, thus neglecting the influence 
of organized spin motion. Although the proliferation of vortices breaks the long range order, 
the spin field can still be organized into domains (where the spins rotate synchronously) 
separated by lines of vortices. The phase transition would be associated in this case, with 
the appearence of a self-organized state, where, although long range order is absent, a kind 
of partial local order, established by separated domains of coherent spins, sets in. 

We consider this possibility and study the spatial distribution of defects around the 
critical temperature. By direct visualization of the spin field rrij = (cos ^j, sin 6'j), we noted 
that the number of isolated defects is negligible, and that they often appear to be in small 
clusters along the domain borders (this may already be seen in Fig. 7 of Ref. The 
problem one encounters with direct visualization is that the spins move too fast to allow the 
observation of coherent structures. To gain some further insight we introduce then a new 
diagnostic to correctly visualize the spatial structures. We update each spin of the lattice 
with one fifth of the sum of itself and its four nearest neighbors and repeat this operation a 
few times, to obtain an effective local magnetization centered on the considered spin. This 
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FIG. 8. Spatial distribution of the local magnetization, showing the domains and the 
interfaces (see text), (a) and (c) orientation (sin^(0^/2)), one passes from white to black 
by a rotation of 7r/4; (b) and (d) intensity (|m|'^|), black are disordered regions, (a) and 
(b): T = 1.02, (c) and (d): T = 1.29. 



local average magnetization is then defined by the iteration: 



m 



(n+l) 



m 



(n) 



+ 



4 

E(n) 



1,...,/, 



where j{i) are the four neighbours of site i, and the number of iterations is typically / = 10. 
The resulting field mf^ is smoother than the initial field mf'^; at these temperatures the 
amplitude of spin motion is large, and it defines a direction at a given site, weighted by the 
orientation of the surrounding spins. 

The effective spin field is shown in Fig. p. The upper plots are linear grey scale images of 
the quantity sin^ (26',f^') (see also ^^^), where Of^ is the angle of the local magnetization, and 



the bottom ones are images of the quantity |nif^|. In the upper plots we can easily locate 
the vortices by looking at pinching of the darker and brighter areas. They appear to be all 
bounded in dipoles or chains of dipoles. The presence of these chains can be interpreted as 
the birth of interfaces separating domains with local order. These latter are more visible 
in the bottom images, where the brighter regions characterize strong local magnetization, 
and therefore locally aligned spins; while the darker regions, where the orientation of spins 
change rapidly with position, indicate the presence of vortex defects and disorder. In the 
left images, representing a system of = 128^ spins at a temperature slightly above the 
critical temperature, T = 1.02, we see that the disordered regions are highly concentrated 
and that they tend to connect themselves along lines. This tendency is even more clear 
in the right images, where we show a system at T = 1.29. We also note that the size of 
the ordered regions, the domains, are much smaller at higher temperatures. At T = 1.02 
the domains occupy connected regions of a size comparable to that of the whole system, at 
higher temperature the domains are instead confined in localized regions. 
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FIG. 9. Typical temporal behavior of a rotator at high temperature. Top, spin 6i(t) 
at the central site for T = 4.100, measured in units of 27r. The spin rotates in the same 
direction and with an almost steady speed for long periods of time. Bottom, plot of 9i(t) 
showing that its fluctuations are fast compared to the evolution of 6i. 



At higher temperature the weight of random regions increases, and the size of domains is 
reduced. In this regime, since the density of potential energy is bounded, most of the energy 
is kinetic and therefore most of the spins are fully rotating, with similar amplitudes. This is 
illustrated in Fig. where the typical temporal behavior of one spin in a high temperature 
field is shown. On the contrary, at low temperatures, angles remain bounded in time. On the 
other hand the spin momentum takes essentially constant values (positive or negative ones, 
depending on the sense of rotation). The individual behavior of the spins also reflects into 
the statistical properties of the system. The angle variance is bounded at low temperatures, 
in agreement with the random phase approximation, while it steadily increases in time at 



high temperatures, as we can see in Fig. |T0|. In analogy with the properties of a perturbed 
pendulum, the change of the topology of trajectories is related to the crossing of a separatrix. 
Below a certain value of the energy density he {he = hiTxr) ~ 1 in our model), the spins 
are collectively trapped and after the transition, as the behavior of the topological defects 
suggests, they start to rotate in organized domains, whose size progressively decreases as 
the energy density increases. 

An important consequence of these results is that the system, at these temperatures, can 
be represented by local ordered regions of synchronized spin motion, and random interfaces 
where the defects accumulate. To these rotating spin domains one may add some level of 
fast (with respect to the time scale of rotation) fluctuations characterized by an effective 
low temperature T, which takes into account the potential energy. In this context we may 
introduce relevant variables Sj, which take the values ±1 according to the sign of spin 
rotation, 6'j. It is important to keep in mind that the system is considered to be in its 
thermodynamical state, and that the temporal average of spin rotation velocity 9f = T is, 
using ergodicity, independent of the spin site. We can then assume that the angle velocity 
is of the form 6i = ili ^ Si^/T, and write 

9i = SiQt + e^it/e) , (4.1) 
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where Q = a/T, and e characterize the fast temporal fluctuations of the function i9j, i.e. 
|e| <C 1 and 16*4^/^)1 ^ SiQt (see Fig. ^). The field 9i is in fact similar to the 9i used in 
but associated to another temperature T: in a rotating reference system, attached to 
each domain (where all the Sj's are equal), it is an effective "low temperature" XY spin field. 
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FIG. 10. The spatial variance of 9 versus time is plotted for different temperatures. 
At the bottom we have the plot for T = 0.4783, while on top it is for T = 2.7905. The 
numerical data is obtained for a lattice size of = 256'^. We clearly see that the variance 
of 9 is an increasing function of time only for T > Tkt- 

The introduction of the new field Si = ±1, which label each lattice site by the sign of the 
spin rotation, suggests an analogy between the high temperature XY system and the Ising 
model. We will exploit this analogy to construct a model aimed at describing the complex 
properties of the system near the Kosterlitz-Thouless transition, which do not resume to the 
simple unbinding mechanism of dipole vortices. 

In order to map the XY system into an Ising system, we take advantage of the time scale 
separation between the relative slow rotation and the fast thermal fiuctuations. Let us now 
average the potential energy over a few rotation periods during which the approximation 
made in ( |4.1|) remains valid. For that purpose, we are reduced to compute the average of 
terms of the form: 

cos(^j — 9j) = cos[(sj — Sj)Qt] cos{9i — 9j) — sin[(sj — Sj)Qt] sm{9i — 9j) , (4.2) 
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where i and j denote two neighbors. The fast temporal behavior of the ^j's assures its fast 
thermahzation at a temperature T and decorrelates the two terms in each product. This 
allows us to split the time average in two steps, first on the fast time dependence and then 
on the slower time scale. The first average leads to 

v{f) 



cos(^i - Oj) = , sin(^i - e,) = 



where we use the results of Section |T|, due to the effective low temperature T within a 
domain. Then, the mean of cos[(sj — Sj)Qt\ is zero if the two spins do not rotate in the 
same direction (do not belong to the same domain) and one if the two spins corotate, which 
means that the average of this terms can be represented by a Heaviside function Q{siSj). 
We then obtain the expression 

^ = ^ E cos (e, - e,) = ^ E ©(^^^.O ^ (4-3) 

<id> <i,j> 

for the mean interaction energy, where the summation is over all neighbors, and the Heaviside 
function just states that only synchronized, corotating spins contribute to the averaged 
potential energy. Using the identity 6(sjSj) = (1 + SjSj)/2, we finally obtain a new effective 
Hamiltonian, 

ff. = iv(I + 2-!fl)-4) E^,..,. (4.4) 

We notice that this effective Hamiltonian is the sum of two different terms, which we denote 
Ht and Hj. The first term Ht depends on the temperature, the size of the lattice and v{T), 
and can therefore be considered as an energy reference term. The second term Hj introduces 
a coupling energy between the spins, and can be recognized as a ferromagnetic Ising-like 
Hamiltonian whose coupling constant Ji{T) = v{T)/4 is a function of the temperature. 
Using an Ising terminology, we have now a new population of spins Sj = ±1, interacting 
with their close neighbors on a square lattice. Therefore, under the assumption of ergodicity 
and taking into account that spins rotate on a well separated time scale with respect to 
their fluctuations, we have mapped the XY system into a Ising-like model; these two models 
are linked through a coupling constant dependent on a temperature related to fluctuations. 
We shall note, that the time dependence of the hamiltonian is masked, but the Sj's are 
still functions of time, and we still are within the microcanonical ensemble. Invoking again, 
the thermodynamical equilibrium, and the fact that both the microcanonial and canonical 
approaches leed to the same thermodynamic limit, we can now continue our study within 
the canonical ensemble. 

As it is well known, the Ising model undergoes a phase transition and generate a spon- 
taneous magnetization Mj when Jj/T > jSjc ~ 0.44 [Q. In order to know if the present 
Ising system reaches the transition, we have to investigate the behavior of v{T). The pre- 
viously observed domains (Fig. |^), are in agreement with the phenomenology of the Ising 
model above its transition, the maximum value of the coupling constant being v(T)/4 < 1/2 
implies that both models are surely in their high temperature states for T > 1/2/3/c ~ 1.14. 
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In the low temperature regime, none of the spins are rotating, T = T and Vj = v, all the 
potential energy is due to the waves. We know that the energy is continuous through the 
Kosterlitz-Thouless transition and also through the Ising transition. As a consequence, we 
expect that v{Txt) = "^iT) at the transition temperature. Moreover, the effective domain 
temperature T has to be an increasing function of T, and then, vj is a decreasing function 
of T, in order to obtain h = T/2 + 2 a.t high temperature, v{t) has to vanish as T ^ cxd. 
The equation 

Ji/T = v{f)/{4T) = (3j, (4.5) 

has therefore a unique solution. A calculation of Tkt, using equations ([4.5| ) and ( p.l8|) gives 
Tkt ~ 0.855, which is in good agreement with the numerical value of Tkt ~ 0.898 found in 
the literature. 

There is however, an important difference between the usual Ising model, and the present 
one, derived from the XY system. It is the global constraint associated to the conservation 
of the total momentum. The mean value of the total momentum P = 0, is given by 
= J^^i = J2 Si- But this last expression is precisely the Ising magnetization: 

^/ = ^E^^ = 0' (4-6) 

i 

which means that only symmetric distribution of positive and negative spins are allowed. 
This constraint prevents the system from undergoing, as in the Ising model, a second order 
phase transition, with spontaneous appearence of a macroscopic magnetization below the 
critical temperature. In the present case, this would mean that a significant fraction of the 



spins would rotate in a preferred direction, but ( |4.6| ) forbids this kind of phase transition, 
and forces the system to accomodate to a vanishing magnetization. In fact what happens 
is that, below the transition temperature, the spins cannot fully rotate, domains dissapear, 
and the long range order is established by a reorganization of spin motions in the form of 
spin waves, which can still generate a nonintensive finite-size magnetization. 



V. CONCLUSION 

In this paper we took a dynamical point of view to analyze the statistical properties of 
the XY model. This approach has the advantage of offering a natural physical framework. 
From the structure of the evolution equations of the spins (which become coupled oscillators 
in the linear approximation), one is prompted to consider the phonons (propagating waves 
in the lattice) as the basic excitation at low temperature. In the thermodynamic limit, and 
assuming that the system reaches an ergodic equilibrium state, it is reasonable to introduce 
random phases in the waves. Using this basic mechanism, we obtained a nonlinear dispersion 
relation containing the fundamental physics of the system. Macroscopic quantities, such as 
the energy and the (finite size) magnetization, are well described by this method up to 
temperatures close to the Kosterlitz-Thouless critical temperature. 

On the other hand, for high temperatures, our dynamical approach allowed us to map 
the XY model into an Ising model with a coupling constant depending on the temperature. 
This relation to an exact model, is beneficial to understand the physics near the transition. 
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and moreover, permits the analytical computation of the critical temperature in excellent 
agreement with Monte Carlo simulations. 

These results lead to a physical description of the Kosterlitz-Thouless transition in terms 
of the change from an ordered state (long range correlation are established by the wave ex- 
citations) to a local ordered state where macroscopic domains are separated by interfaces 
populated with topological defects. Topological defects have then a tendency to bind them- 
selves in clusters. The generation of disorder by the unbinding of dipoles, is accompanied 
by the fragmentation of the low energy unique domain of non rotating spins into separated 
regions of synchronized rotating spins. 
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